Batch Calibration Thermogram

Load packages

library("Thermimage")
library("fields")
library("imager")
library("magrittr")
library("spatstat")

Set configurations

## Noise filering
NF = TRUE
## Histgoram Equalization
HE = TRUE

## INPUT_PATH_FOLDER
input = "../Data/"
## OUTPUT_PATH_FOLDER
output = "4_HE_NF/"

## Gassian sigma value for noise filtering
sigma = 0.8

## Situational parameters
OD = 10
RH = 53
AtmosT = 28

Load path images from directory

filenames <- list.files(input, pattern="*.JPG", full.names=FALSE)
filenames

  [1] "DJI_0001_R.JPG" "DJI_0003_R.JPG" "DJI_0005_R.JPG" "DJI_0007_R.JPG" "DJI_0009_R.JPG" "DJI_0011_R.JPG"
  [7] "DJI_0013_R.JPG" "DJI_0015_R.JPG" "DJI_0017_R.JPG" "DJI_0019_R.JPG" "DJI_0021_R.JPG" "DJI_0023_R.JPG"
 [13] "DJI_0025_R.JPG" "DJI_0027_R.JPG" "DJI_0029_R.JPG" "DJI_0031_R.JPG" "DJI_0033_R.JPG" "DJI_0035_R.JPG"
 [19] "DJI_0037_R.JPG" "DJI_0039_R.JPG" "DJI_0041_R.JPG" "DJI_0043_R.JPG" "DJI_0045_R.JPG" "DJI_0047_R.JPG"
 [25] "DJI_0049_R.JPG" "DJI_0051_R.JPG" "DJI_0053_R.JPG" "DJI_0055_R.JPG" "DJI_0057_R.JPG" "DJI_0059_R.JPG"
 [31] "DJI_0061_R.JPG" "DJI_0063_R.JPG" "DJI_0065_R.JPG" "DJI_0067_R.JPG" "DJI_0069_R.JPG" "DJI_0071_R.JPG"
 [37] "DJI_0073_R.JPG" "DJI_0075_R.JPG" "DJI_0077_R.JPG" "DJI_0079_R.JPG" "DJI_0081_R.JPG" "DJI_0083_R.JPG"
 [43] "DJI_0085_R.JPG" "DJI_0087_R.JPG" "DJI_0089_R.JPG" "DJI_0091_R.JPG" "DJI_0093_R.JPG" "DJI_0095_R.JPG"
 [49] "DJI_0097_R.JPG" "DJI_0099_R.JPG" "DJI_0101_R.JPG" "DJI_0103_R.JPG" "DJI_0105_R.JPG" "DJI_0107_R.JPG"
 [55] "DJI_0109_R.JPG" "DJI_0111_R.JPG" "DJI_0113_R.JPG" "DJI_0115_R.JPG" "DJI_0117_R.JPG" "DJI_0119_R.JPG"
 [61] "DJI_0121_R.JPG" "DJI_0123_R.JPG" "DJI_0125_R.JPG" "DJI_0127_R.JPG" "DJI_0129_R.JPG" "DJI_0131_R.JPG"
 [67] "DJI_0133_R.JPG" "DJI_0135_R.JPG" "DJI_0137_R.JPG" "DJI_0139_R.JPG" "DJI_0141_R.JPG" "DJI_0143_R.JPG"
 [73] "DJI_0145_R.JPG" "DJI_0147_R.JPG" "DJI_0149_R.JPG" "DJI_0151_R.JPG" "DJI_0153_R.JPG" "DJI_0155_R.JPG"
 [79] "DJI_0157_R.JPG" "DJI_0159_R.JPG" "DJI_0161_R.JPG" "DJI_0163_R.JPG" "DJI_0165_R.JPG" "DJI_0167_R.JPG"
 [85] "DJI_0169_R.JPG" "DJI_0171_R.JPG" "DJI_0173_R.JPG" "DJI_0175_R.JPG" "DJI_0177_R.JPG" "DJI_0179_R.JPG"
 [91] "DJI_0181_R.JPG" "DJI_0183_R.JPG" "DJI_0185_R.JPG" "DJI_0187_R.JPG" "DJI_0189_R.JPG" "DJI_0191_R.JPG"
 [97] "DJI_0193_R.JPG" "DJI_0195_R.JPG" "DJI_0197_R.JPG" "DJI_0199_R.JPG" "DJI_0201_R.JPG" "DJI_0203_R.JPG"
[103] "DJI_0205_R.JPG" "DJI_0207_R.JPG" "DJI_0209_R.JPG" "DJI_0211_R.JPG" "DJI_0213_R.JPG" "DJI_0215_R.JPG"
[109] "DJI_0217_R.JPG" "DJI_0219_R.JPG" "DJI_0221_R.JPG" "DJI_0223_R.JPG" "DJI_0225_R.JPG" "DJI_0227_R.JPG"
[115] "DJI_0229_R.JPG" "DJI_0231_R.JPG" "DJI_0233_R.JPG" "DJI_0235_R.JPG" "DJI_0237_R.JPG" "DJI_0239_R.JPG"
[121] "DJI_0241_R.JPG" "DJI_0243_R.JPG" "DJI_0245_R.JPG" "DJI_0247_R.JPG" "DJI_0249_R.JPG" "DJI_0251_R.JPG"
[127] "DJI_0253_R.JPG" "DJI_0255_R.JPG" "DJI_0257_R.JPG" "DJI_0259_R.JPG" "DJI_0261_R.JPG" "DJI_0263_R.JPG"
[133] "DJI_0265_R.JPG" "DJI_0267_R.JPG" "DJI_0269_R.JPG" "DJI_0271_R.JPG" "DJI_0273_R.JPG" "DJI_0275_R.JPG"
[139] "DJI_0277_R.JPG" "DJI_0279_R.JPG" "DJI_0281_R.JPG" "DJI_0283_R.JPG" "DJI_0285_R.JPG" "DJI_0287_R.JPG"
[145] "DJI_0289_R.JPG" "DJI_0291_R.JPG" "DJI_0293_R.JPG" "DJI_0295_R.JPG" "DJI_0297_R.JPG" "DJI_0299_R.JPG"
[151] "DJI_0301_R.JPG" "DJI_0303_R.JPG" "DJI_0305_R.JPG" "DJI_0307_R.JPG" "DJI_0309_R.JPG" "DJI_0311_R.JPG"
[157] "DJI_0313_R.JPG"

Read image function

# Create a function with arguments.
new.read_image <- function(image_file, 
                           input_path="", output_path="Output/",
                           NF=FALSE, 
                           sigma=0.8,
                           HE=FALSE,
                           OD=NULL,
                           RH=NULL,
                           AtmosT=NULL) {
  ## Create image path
  image_path = paste(input_path, image_file, sep="")

  ## Extract meta-tags from thermal image file ##
  cams<-flirsettings(image_path, exiftool="installed", camvals="")

  ## Set variables for calculation of temperature values from raw A/D sensor data
  Emissivity<-cams$Info$Emissivity      # Image Saved Emissivity - should be ~0.95 or 0.96
  ObjectEmissivity<-0.96                # Object Emissivity - should be ~0.95 or 0.96
  dateOriginal<-cams$Dates$DateTimeOriginal
  dateModif<-   cams$Dates$FileModificationDateTime
  PlanckR1<-    cams$Info$PlanckR1                      # Planck R1 constant for camera
  PlanckB<-     cams$Info$PlanckB                       # Planck B constant for camera
  PlanckF<-     cams$Info$PlanckF                       # Planck F constant for camera
  PlanckO<-     cams$Info$PlanckO                       # Planck O constant for camera
  PlanckR2<-    cams$Info$PlanckR2                      # Planck R2 constant for camera
  ATA1<-        cams$Info$AtmosphericTransAlpha1        # Atmospheric attenuation constant
  ATA2<-        cams$Info$AtmosphericTransAlpha2        # Atmospheric attenuation constant
  ATB1<-        cams$Info$AtmosphericTransBeta1         # Atmospheric attenuation constant
  ATB2<-        cams$Info$AtmosphericTransBeta2         # Atmospheric attenuation constant
  ATX<-         cams$Info$AtmosphericTransX             # Atmospheric attenuation constant
  if (is.null(OD)) {
    OD<-          cams$Info$ObjectDistance                # object distance in metres
  }
  FD<-          cams$Info$FocusDistance                 # focus distance in metres
  ReflT<-       cams$Info$ReflectedApparentTemperature  # Reflected apparent temperature
  if (is.null(AtmosT)) {
    AtmosT<-      cams$Info$AtmosphericTemperature        # Atmospheric temperature
  }
  IRWinT<-      cams$Info$IRWindowTemperature           # IR Window Temperature
  IRWinTran<-   cams$Info$IRWindowTransmission          # IR Window transparency
  if (is.null(RH)) {
    RH<-          cams$Info$RelativeHumidity              # Relative Humidity
  }
  h<-           cams$Info$RawThermalImageHeight         # sensor height (i.e. image height)
  w<-           cams$Info$RawThermalImageWidth          # sensor width (i.e. image width)

  ## Rotate image
  temp_image = rotate270.matrix(readflirJPG(image_path))

  ## Convert raw to temp data
  temp_image = raw2temp(temp_image, ObjectEmissivity, OD, ReflT, AtmosT, IRWinT, IRWinTran, RH,
                        PlanckR1, PlanckB, PlanckF, PlanckO, PlanckR2)
  
  ## Noise filtering with gauss filter
  if (NF == TRUE) {
    temp_image =  as.im(temp_image)
    temp_image = blur(temp_image, sigma, bleed=FALSE)
    temp_image =  as.matrix(temp_image)
  }

  
  # hist(temp_image, main=paste("Histogram ", image_file))
  ## Save temperature image
  saveRDS(temp_image, file = paste(output_path, tools::file_path_sans_ext(image_file), ".RData", sep=""))

  ## Return either min or max or every temperature value
  if (HE == TRUE) {
    temp_image
  } else {
    min <- min(temp_image)
    max <- max(temp_image)
    values <- c(min, max)
  }
}

Load images

glob_temp_range <- unlist(lapply(filenames, new.read_image, 
                                 input_path=input,
                                 output_path=output,
                                 NF=NF,
                                 sigma=sigma,
                                 HE=HE,
                                 OD=OD,
                                 RH=RH,
                                 AtmosT=AtmosT))

Calculate global min and max values

min_temp <- min(glob_temp_range)
max_temp <- max(glob_temp_range)

Calculate ECDF

if (HE == TRUE) {
  f <- ecdf(glob_temp_range)
  hist(glob_temp_range, main=paste("Global Histogram"))
  hist(f(glob_temp_range), main=paste("Global Equalized Histogram"))
} else {
  f <- NULL
}

Normalize images

new.normalize <- function(image_file, output_path="4_HE_NF/",
                          min, max, 
                          HE=FALSE, f_ecdf=NULL) {
  ## Create image path
  image_path = paste(output_path, tools::file_path_sans_ext(image_file), ".RData", sep="")
  temp_image <- readRDS(image_path)
  # unlink(image_path)
  
  ## Normalize data
  temp_image <- mirror.matrix(temp_image)

  if (HE == TRUE) {
    equalized <- f_ecdf(temp_image)
    equalized <- as.cimg(equalized, dim=dim(temp_image)) 
    plot(equalized, rescale=FALSE, main=paste("Histogram equalized", image_file))
    save.image(equalized, file=paste(output_path, image_file, sep=""), quality = 1)
  } else {
    temp_image <- (temp_image - min) / (max - min)
    # hist(x=temp_image, main=paste("Histogram of normalized values", image_path), xlab='normalized value')

    ## Plot and save temperature image
    temp_image = as.cimg(temp_image, dim=dim(temp_image))
    save.image(temp_image, file=paste(output_path, image_file, sep=""), quality = 1)
    plot(temp_image, rescale=FALSE)
  }
}

Plot normalized images

temp = lapply(filenames, new.normalize, 
              output_path=output, 
              min=min_temp, max=max_temp, 
              HE=HE, 
              f_ecdf=f)